function f=vectorField_CRTBP(y,G,mu)

%the distances
r1=sqrt((mu+y(1))^2+(y(2))^2+(y(3))^2);
r2=sqrt((y(1)-(1-mu))^2+(y(2))^2+(y(3))^2);
%masses
m1=1-mu;
m2=mu;

f=[y(4); 
    y(5); 
    y(6); 
    y(1)+2*y(5)+G*m1*(-mu-y(1))/(r1^3)+G*m2*(1-mu-y(1))/(r2^3); 
    y(2)-2*y(4)-G*m1*(y(2))/(r1^3)-G*m2*y(2)/(r2^3); 
    -G*m1*y(3)/(r1^3)-G*m2*y(3)/(r2^3)];